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ABSTRACT 



We compare various approaches to find the most efficient method for the practical computation of the lightcurves (integrated bright- 
nesses) of irregularly shaped bodies such as asteroids at arbitrary viewing and illumination geometries. For convex models, this 
reduces to the problem of the numerical computation of an integral over a simply defined part of the unit sphere. We introduce a 
fast method, based on Lebedev quadratures, which is optimal for both lightcurve simulation and inversion in the sense that it is the 
simplest and fastest widely applicable procedure for accuracy levels corresponding to typical data noise. The method requires no 
tessellation of the surface into a polyhedral approximation. At the accuracy level of 0.01 mag, it is up to an order of magnitude faster 
than polyhedral sums that are usually applied to this problem, and even faster at higher accuracies. This approach can also be used 
in other similar cases that can be modelled on the unit sphere. The method is easily implemented in lightcurve inversion by a simple 
alteration of the standard algorithm/software. 

Key words. Methods: numerical - Techniques: photometric - Scattering - Minor planets, asteroids: general 



1. Introduction 

Many quantities in astronomy, physics, and other sciences are 
most naturally represented as functions on the unit sphere S 2 . 
Typical examples in astrophysics are the celestial sphere or the 
surfaces of celestial bodies. A standard operation for such func- 
tions is to compute their integrals over the whole sphere or parts 
of it, usually with some weight functions, to express total sums, 
averages, moments, and so on. Integrals over any shape other 
than the sphere, but ones that can be mapped one-to-one onto 
S 2 , are handled using the surface element (Jacobian) as a weight 
function and are thus reduced to exactly the same spherical prob- 
lem. 

Brightness integrals or lightcurves are an important example 
of these surface integrals. As shown in Durech & Kaasalainen 
(2003 ), virtually all lightcurves of single asteroids in realistic 
observing geometries can be explained, i.e., fitted down to noise 
level, by convex surfaces approximating the global shape of the 
target. Thus, for all photometric purposes, such targets can be 
represented as convex models. Convex models require no ray- 
tracing as their brightnesses can be expressed as integrals over a 
simply defined part of S 2 . What is more, most models are very 
conveniently and compactly described by function series (based 
on spherical harmonics) in lightcurve inversion (Kaasalainen et 
al.[l992a 2001 2006). One can thus, in principle, compute their 
total brightness (especially that of models of real asteroids) at ar- 
bitrary viewing geometries without the polyhedral surface repre- 
sentations that are the standard means of integration in lightcurve 
simulation and inversion. In this paper, we show that function 
integration is faster and, in many ways, simpler than polyhedral 
modelling. This is particularly advantageous for large-scale sim- 
ulations and inverse problems involving thousands of targets and 
the repeated computation of millions of lightcurve points. 



A traditional way to compute surface integrals is to tessellate 
the surface, i.e., to represent it as a polyhedral approximation, 
usually with triangular facets of similar size (triangulation). The 
integral is then approximated as a sum over the facet areas multi- 
plied by the corresponding values of the integrand functions. The 
approximation can be made arbitrarily accurate by making the 
triangulation mesh denser. However, in the (log N, log A)-plane, 
where A is the approximation error and N the number of facets, 
i.e., of function evaluations required, the accuracy log A(logA0 
improves only with a linear slope of - 1 . To reduce the error by 
a half, one must double the number of computations. This is 
the least efficient way of computing the integrals, equivalent to 
the elementary way of computing one-dimensional integrals as 
Riemann sum approximations by dividing the integration inter- 
val into equally large steps. 

In the one-dimensional case of real-valued functions on R, or 
its Cartesian multidimensional version R", the standard way of 
improving the computation is to use quadratures such as Gauss- 
Legendre (or Gauss-Chebyshev etc., depending on weight func- 
tions). For large n, the only computationally feasible way is 
to use Monte Carlo techniques. The sphere S 2 has a markedly 
different topology from R 2 ; using Gaussian quadratures, some 
of the computations are wasted on integration points clustered 
too tightly around the two poles, and the integration scheme is 
thus manifestly dependent on the coordinate frame. Evidently, 
the most natural quadrature for S 2 is quite different from the 
Cartesian case. 

Lebedev quadrature (Lebedev & Laikov 1999) is designed 
to be rather automatically optimal for smooth and continuous 
integrand functions over the whole of S 2 , in the same sense that 
Gaussian quadrature is optimal for R. The Lebedev scheme com- 
bines the almost equidistant and rotationally invariant distribu- 
tion of points on S 2 with the efficiency of quadratures (a low ver- 
tical positioning and a slope steeper than -1 for the log A(log N)- 
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curve). However, it is not designed for integrals only over a part 
of S 2 . For these, we need to consider and test various approaches 
to see which is the most useful one in practice. Our purpose is to 
find the optimal approach for accuracy levels corresponding to 
typical data noise. 

This paper is organized as follows. In Section 2 we review 
the problem of the computing of lightcurve or brightness inte- 
grals on S 2 , and discuss its numerical representation in the form 
of quadratures. In Section 3 we present efficiency results from 
a various sample cases, and also note a particular formulation 
for integrals on triaxial ellipsoids. In Section 4 we describe the 
use and efficiency of the quadrature scheme in both inversion 
and simulations (e.g., Monte Carlo sampling). In the concluding 
Section 5 we sum up and discuss some aspects of the quadrature 
approach. 

2. Brightness integrals as quadratures 

A detailed discussion of the theoretical aspects of the problem 
of computing the total brightness of a convex surface is given 
in Kaasalainen et al. (1992a, 20011 120061 ). so here we only re- 
view the main points and notations relevant to our search of the 
optimal method. 

Let the illumination (subsolar) and viewing (sub-Earth) di- 
rections in a coordinate system fixed to the target body be de- 
noted by a>o, oj e S 2 . Here entities in S 2 , defined by two direc- 
tion angles, are identified with unit vectors in R 3 . Thus, e.g., the 
outward unit normal vector 77 e S 2 is given by r\ — rj(§, iff) (with 
# measured from the pole); i.e., 



rji = sin)? cos iff, rj2 = sin)? sin iff, 773 = cos#. 



(1) 



Denoting the portion of the surface both visible and illuminated 
by A + , the total (integrated) brightness of the body is given by 



L((i)o,to)= I S(fj.,fio,(x)G(d-,i[f)dcr, 
Ja + 



where 

dcr = sin & d§ dift, 



(2) 



(3) 



and the scattering function SQ-i,fio,a) is taken to depend on the 
viewing geometry quantities 



fx = oj • 77, fig = a>o ■ t], cos a - u)q ■ u). 



(4) 



Thus A + includes those points on the surface for which /iq, ix > 0. 
The curvature function G > of the surface is 



G(#, iff) = 



/()?, if/) 



sin)? 

where J := \J\ is the norm of the Jacobian vector 
w ,„ „ dx dx 



(5) 



where x(tf, if/) gives the surface as a function of the surface nor- 
mal direction. 

The standard choice, adopted in Kaasalainen et al. (120011) . 
for the numerical computation of (f2]i is the tessellation of S 2 
into N almost equal-sized facets approximating the surface as a 
polyhedron: 



L(« , oj) * J] S A^. «) <KVj) Acr„ 



where rjj is the surface normal of the facet j, and the tessellation 
is done on S 2 with any standard scheme, with Act, as the area of 
the facet j on S 2 . Thus the area of the facet j on the actual sur- 
face is G(rij) Actj, and this is solved for in the inverse problem, 
leading to surface reconstruction via the Minkowski procedure 
(Kaasalainen et al. I1992al [200TT>. 

Another, more analytical option for evaluating (f2| is to use 
quadratures rather than the geometrically intuitive equal-facet 
tessellation. The main principle in any quadrature scheme is to 
choose the evaluation points of the integrand function (and their 
corresponding weights) in a cleverer way than a brute-force ap- 
proximation of a Riemann sum. This choice is based on an as- 
sumption of the form of the function. Gaussian quadrature is 
based on the fact that any one-dimensional integral (with finite 
itegration limits) over any polynomial fix) of (2N - 1 )th degree 
can be evaluated exactly by choosing a set, independent of / and 
tabulated in one step, of merely TV unevenly distributed evalua- 
tion points or abscissae jc, (to be scaled linearly according to the 
integration interval), and their weights w, so that 



J fix) dx = ^ f(Xi)Wi 



(8) 



(Arfken 1985), so the sum can be used as a good approxima- 
tion of integrals over any other well-behaved functions fix) (de- 
pending on the number of abscissae chosen). Consequently, any 
two-dimensional integral 



a Jc(v) 



fiu, v) du dv 



(9) 



over a double polynomial fiu, v) of degrees L, M in u, v can be 
evaluated by using (L + 1)(M + l)/4 tabulated abscissae and 
their weights. Here we are interested in the spherical harmon- 
ics F™(m, v) that can be expressed as double polynomials in u, v. 
Since the v-part is expressed by e" m \ i.e., two polynomials from 
the functions sin mv and cos mv, and we usually have M = L, 
the maximum degree of the spherical harmonics used, the num- 
ber of abscissae needed for a Gaussian quadrature over spherical 
harmonics is (L + l) 2 /2. The quadrature is usually a combina- 
tion of Gauss-Legendre and Gauss-Chebyshev to account for the 
trigonometric functions with the corresponding weight function. 

We now investigate the available options for writing Eq. (f2j) 
in a standard quadrature form. In Kaasalainen et al. (11992at . the 
analytical handling of Eq. (f2]), leading to the uniqueness proof 
for the inverse problem, required transforming the integrand into 
a coordinate system in which ojq, oj are in the new xy-plane. This 
yields simple integration limits, but comes at the price of a more 
complicated integrand: 



L( oj q , oj) 



where 



Pr 

a JO 



ioJo, w)G()?, i//)S (/J, /io, a) dcr, 



(6) yD = sin iff sin §, fio = sin(</> - a) sin §, 



(10) 



(ID 



and the operator P«(wo, oj) transforms a function into the new 
system obtained by the frame rotation matrix R(wo, oj)- Thus, 

P r (oj , oj)G(&, iff) = G(FT Vo, w)iK*. <A)X (12) 

or, if G is of the form G = Gi{Yf l i§, iff)}), where {Y™ (■&, iff)} de- 
notes some set of spherical harmonics with a number of various 
degrees and orders /, m, 



(7) P R (oj ,oj)Gi&,iff) = G 



J]r' i m ioj ,oj)Yii&,iff)\ 



(13) 
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where r'™(ajo, &>) are me so-called D-matrices of spherical har- 
monics (Kaasalainen et al. I1992al l. 

The integration limits of Eq. ( TT0T > nominally allow two- 
dimensional Gaussian quadrature, but the forms (IT2b and (TT~3T > 
both cause significant computational overhead, as well as com- 
plications in inverse problem applications. Using Eq. (TT~2t calls 
for additional computations as the evaluation points of G(rj) 
are redefined, and the transformation (13[ introduces additional 
multiplications for function values at the evaluation points. 
Furthermore, the inverse problem solution typically requires the 
derivatives of L with respect to the rotation parameters of the 
target that define ojq, u>, and this leads to complicated derivative 
chains in both fl% and dT3l > which increase the computational 
cost further. Note also that, due to the existence of S in the in- 
tegrand, and the fact that we usually use functions exp(y ; m ) to 
guarantee positivity, the nominal number of quadrature points 
discussed above is not sufficient for a series of degree L. 

Still another way to do Gaussian quadrature over a part of S 2 
is to keep the integrand simple but to write the integration limits 
as, e.g., 



Distribution of Lebedev points(Vertices) (N=230) 



L= I I S(ii,fi ,a)G(-»,i//) sm$d§dif/, (14) 



but now the boundary equations fi = 0, fio = are implicit for 
both if/, and ftiiifr), so setting the Gaussian grid calls for compu- 
tationally expensive root-finding, and the parameter derivatives 
are also somewhat cumbersome, as above. 

Further overhead is caused in typical repeated lightcurve 
computations when the function G(rf) cannot be evaluated at a 
set of fixed 77,- in one step. This happens in the above integrals 
due to the changing viewing geometry that alters the integration 
limits, hence the sampling points of G(rj) for each computation 
of L; for Eq. (TlOb this is caused by a(a>o, u>) in the integration 
limits. If the number of lightcurve points for one target is M, 
that of computations required for evaluating S (fi, fio, ct) is I, and 
K function evaluations are needed for determining G(tj) (with, 
e.g., spherical harmonics), then the computational cost of fixed- 
G(rji) evaluation grows as NMI, while that of Gaussian quadra- 
ture grows as NM(I + K), when N is the size of the 77,-grid. Since 
typically K » /, the K/I-fold overhead is considerable. 

In view of the above, we thus discard the use of the above 
forms for rendering (f2]i to a form suitable for Gaussian quadra- 
ture as the potentially reduced number of evaluation points, as 
compared to tessellation, is not nearly sufficient to compensate 
for the computational overhead in the applications in practice 
we have in mind. In these, the accuracy levels are around 1% 
(0.01 mag), where a few hundred facets suffice for tessellation 
computations in any case. In applications where high accuracy 
is needed, Gaussian quadrature over a part of S 2 via Eqns. ([Tol l 
or (fl4l can be useful as the efficiency of any quadrature grows 
very fast with growing N when the integrand is smooth and con- 
tinuous over the integration region. 

Quadratures can also be introduced by formally writing the 
integral (f2]) over the whole of S 2 by introducing a transformation 
zeroeing S(fi,/uo, a) in places: 



, 0, n < or no < 

*W,[io,a) - < s(jl flo a ^ elsewhere, 



so that we have 



I G(&,il/)S(fi,/jQ,a)dcr. 
Jo 



(15) 



(16) 



0.5, 



-0.5- 




-1 -1 



Fig. 1. The distribution of Lebedev points on the unit sphere. 



Now the integration limits are simple constants, so this can be 
readily computed with quadratures, and the evaluation points are 
the same for each viewing geometry. The only hitch is that the 
integrand is usually no longer smooth or not necessarily even 
continuous everywhere in the integration region (in contrast with 
polynomials), so it is not possible to give a simple estimate of the 
minimum size of the quadrature grid. 

While we could use Gaussian quadrature to evaluate Eq. 
( TTSI l. Lebedev quadrature is a more natural choice. It is not based 
on the double-integral form, but instead evaluates integrals on 
the whole of S 2 ; i.e., 



Jo 



f(u,v) sinudu dv 



as an approximation 

N 

1 ~ ^f (lli > v ') w ''' 



(17) 



(18) 



where the evaluation points (m,, v,) and their weights w, can be 
tabulated for various N with standard algorithms (Lebedev & 
Laikov [T999l l. These tables can be called in the evaluations in 
the same way as with other quadratures. Software for creating 
the tables is available on the Internet Q, and sample tables and 
Matlab software are available from the authors. 

For Lebedev quadrature, the number of required points in 
the grid is 2/3 of the Gaussian case, scaling approximately as 
(L + l) 2 /3 when the integrand / is assumed to be a spherical 
harmonics series with maximum degree L so that the quadra- 
ture yields the exact value of the integral (Lebedev & Laikov 
1999, Wang et al. 2003). Lebedev quadrature is also likely to 
give more accurate values than a Gaussian one in the case of 
a non-smooth integrand because of its almost even distribution 
of points on S 2 that also makes the computation virtually inde- 
pendent of the coordinate frame. This is depicted in Fig.[T] ver- 
tices are the Lebedev points, and the connecting edges are added 
by Delaunay triangulation to illustrate how the node distribution 
somewhat resembles that of typical tessellation. 



1 see, e.g., http://people.sc.fsu.edu/~jburkardt/ 
f_src/sphere_lebedev_rule/sphere_lebedev_rule .html; 
http : //www .mathworks . com/matlabcentral/f ileexchange/ 
27097-getlebedevsphere 



3 



M. Kaasalainen et al.: Optimal brightness computation 



10 J 



Error Curves: Triangulationf) Lebedev(o) 



10 ' 



r *** 



Oo 

o 



10 10 10 10 

Number of Evaluation points 



10" 



10 



Error Curves: Triangulation(*) Lebedev(o) 
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Fig. 2. The error curves log A(log N) for computing the surface Fig. 3. The error curves log A(log AO for computing the area of 
integral ( l20b . Asterisks denote integration by triangulation, while an ellipsoid. Asterisks denote integration by triangulation, while 



circles are for integration by quadratures. 



circles are for integration by quadratures. 



The derivation of the abscissae and weights is more com- 
plicated than in the Gaussian one-dimensional case due to the 
wholesale operation over S 2 that takes rotational symmetries 
into account. However, the tabulated end result is just as easy 
to use as in the one-dimensional Gaussian quadrature. Lebedev 
quadrature is essentially an advanced method for arranging eval- 
uation points such that their distribution is approximately but 
not completely even, and consequently their weights w-, are not 
uniform (in contrast with standard tessellation that aims at sim- 
ilar Act). This predictive use of symmetry gives a competitive 
edge over both equal-facet tessellation and Gaussian quadrature. 
Most of the points have weights of roughly similar sizes, but the 
values of some weights are much smaller than the average and, 
depending on the chosen number of points, some weight values 
may even be negative. 

In fact, for Lebedev schemes with positive weight values, 
Lebedev quadrature is formally equivalent to a customized tes- 
sellation scheme, the weights being interchangeable with the 
facet areas of a polyhedron. Here 2; w i = 4n; by the construc- 
tion of Lebedev quadrature, the polyhedron approximates the 
unit sphere. Even though no exact description of such an equiv- 
alence relation is needed for the problem here, we note that the 
locations of the vertices corresponding to the areas can be de- 
termined by, e.g., solving the polyhedral Minkowski problem 
(Kaasalainen et al. 2001] [2006l >. Figure [1] portrays a dual im- 
age of such a polyhedral tessellation, vertices corresponding to 
facets. Lebedev points are thus not directly the vertices of an op- 
timal tessellation, but they indirectly define one for certain cases. 

Now, with 

ftn) - G{Tj)S , (19) 

many of the f(rji) vanish, so we take the quadrature sum only 
over the points for which hq,h > 0, just as in the tessellation 
sum. 

To illustrate the remarkable power of Lebedev quadrature on 
the whole of S 2 for a smooth and continuous integrand, we con- 
sider the simple integral 

r*2n r*7T ^ 

1=1 I sin ududv = 27iTn3, (20) 

Jo Jo 2 - cos u 



and plot the curves log A(log AO of the relative error A of the nu- 
merical value (from comparison with the analytical one) in the 
(log N, log A)-plane in Fig. [2] The error of the triangulation sum 
is plotted with asterisks (with N for the number of facets), and 
that of the quadrature with circles; the saturation of the quadra- 
ture is due to machine accuracy. This curve also reflects the fact 
that quadratures over a well-behaved integrand always eventu- 
ally win over equal-facet or standard tessellation (here we refer 
to this triangulation by the terms "tessellation" and "triangula- 
tion") even if there is computational overhead, as discussed with 
Gaussian quadrature over a part of S 2 above. 

We are now left with two methods to compare: the tessella- 
tion sum (01 and the Lebedev sum (1181 1. Our hypothesis of the 
superior efficiency of using § with Lebedev quadrature on S 2 
rests on the assumption that the quadrature can handle the non- 
smooth behaviour of S without invoking too many function eval- 
uations; i.e., if we think in terms of an optimal tessellation via 
the Lebedev procedure, it should retain its superior properties 
even when shadowing is included. In the next section we show 
that this is indeed the case. Apparently the fact that most realistic 
scattering models contain the common factor fifio, thus making 
the scattering function S (jd, yUo) vanish at the integration bound- 
ary, also helps to retain the accuracy in the quadrature since the 
integrand is then at least continuous on all of S 2 . 

3. Benchmark examples 

In this section, we present some benchmark tests of the accuracy 
of the methods of computing surface integrals as the size of the 
evaluation grids grows. In all figures, we have denoted the nom- 
inal number (AO of function evaluations in the tessellation and 
quadrature sums. In brightness examples, more than half of the 
evaluation points are skipped in the computation due to // < or 
/A) ^ (i.e., S = 0), but the relative portions of the skipped terms 
are the same for both sums, so the relative difference between the 
numbers of computations stays the same. 

As another example of Lebedev quadrature with smooth in- 
tegrands over the whole of S 2 , we plot in Fig. [3]the relative ap- 
proximation error A of the area A of an ellipsoid against the size 
N of the evaluation grid (asterisks for tessellation facets and cir- 
cles for Lebedev points; A is not analytically computable, but 
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Error Curves with Ellipsoid a=8,b=6,c=5 for 500 tests 
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Error Curves: Tessellation^) Lebedev(o) for 500 tests 
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Fig. 4. Average logA(logA0 for geometrically scattering ellip- 
soids at various random observing geometries (asterisks for tes- 
sellation, circles for quadrature). N is the nominal number of 
evaluation points in a sum (over half of these are skipped). 
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Fig. 5. Average log A(log N) for general asteroid shapes and 
Lommel-Seeliger/Lambert scattering at various random observ- 
ing geometries. N is the nominal number of evaluation points in 
a sum (over half of these are skipped). 



we use a numerical elliptic integral routine for establishing the 
accuracy). We thus compute 



r-ln f7T 

Jo Jo 



G(#,ifr)sm#d#difr, 



(21) 



where G(ft, if/) for an ellipsoid with semiaxes a, b, c is given in 
Kaasalainen et al. (19923: 



G(&,iff) 



abc 



y {a sin -& cos if/) 2 + (b sin ft sin if/) 2 + (c cos ff)' 



.(22) 



Again, the quadrature is clearly superior to tessellation, as ex- 
pected. The saturation at high accuracy is mostly due to numeri- 
cal roundoff effects. 

The situation changes visibly when we consider the bright- 
ness of a geometrically scattering ellipsoid, i.e., S = fi. This can 
be computed analytically (Ostro & Connelly |1984| l for compari- 
son: 



71 

L — —abc 
2 



a) + 



w r Mw 



(23) 



where 
M 



( 1/a 2 (T 
1/b 2 
1/c 2 , 



(24) 



We chose random realistic observing geometries with roughly 
the same proportions of terms skipped in the quadrature (facets 
omitted in the tessellation sum). Now the integrand function S 
is not even continuous over S 2 , and this shows in the quadrature 
error curve in Fig. [4] that is similar to the tessellation curve in 
its slope. However, it lies clearly below the latter, so even now 
its performance is noticeably better for our purposes. By chang- 
ing the viewpoint 90 degrees and looking at the A(A0-curves as 
plots of the inverse function N(A) of the number of evaluations 
required for a given accuracy (especially around the relative er- 
ror 0.01), one can clearly see the difference in computational cost 
between the two approaches. 



Realistic scattering functions make § continuous, which 
turns out to be of an important value. As a more general test, 
we examine the accuracy performance of Lebedev and tessella- 
tion grids by computing the brightness L at a number of random 
observing geometries for general convex shapes (asteroid mod- 
els obtained from real lightcurve data). Now G($, iff) is given 
by exponential spherical harmonics series of shape models as 
used in inversion (Kaasalainen et al. 2001 ). The scattering model 
used here was of the form combining the Lommel-Seeliger and 
Lambert laws, 



S(po,fi,a) = f(a)fi fi 



1 



+ C 



(25) 



which is typically used in lightcurve inversion; i.e., f(a) has no 
influence on the plot in Fig. [5] The reference value of L for com- 
puting A was obtained by applying a dense tessellation grid. 
Even though the performance of the quadrature is not as good 
as for actual integrals over the whole of S 2 , the Lebedev -based 
scheme is clearly more efficient than tessellation also in that its 
error curve has a steeper slope. The same behaviour applies to, 
e.g., ellipsoids. Also, in Fig. [5] there are no scattered points such 
as the few seen in the Lebedev curve in Fig. [4] This is another 
beneficial effect of the continuity provided by the product /j/jq: 
when Lebedev points are sparse, a discontinuous integrand is 
more vulnerable to the geometry of the integration depending 
on whether individual Lebedev points (and their weights) are in- 
cluded in the quadrature sum or not. In general, predicting the 
performance of the quadrature-based approximation for discon- 
tinuous integrands, especially with a low number of evaluation 
points, is usually possible only by numerical tests. 

The choice of the surface parametrization determines how 
the quadrature evaluation points are distributed, and for a priori 
known convex shapes we can, of course, use any representations 
of the surface points x(u, v), not necessarily just the Gaussian 
image x(rf). In such cases we write 



L{u) , a/) 



-j 



S (fio, fi, a) J(u, v) du dv, 



(26) 
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and the normal vector required for determining /j,/jo is 77 = J /J. 
For example, ellipsoidal surfaces can be given in the standard 
parametrization 



x\ — a sin u cos v, X2 — b sin u sin v, X3=ccosm. 



(27) 



Thus, for an ellipsoid, we have, using u, v instead of §, iff in (O, 

— ' - = yj(bc sin u cos v) 2 + (ac sin u sin v) 2 + (a/> cos m) 2 . (28) 
sin m 

Since < u < n and < v < 2n, the integral (|26*T i is reduced 
to the standard Lebedev quadrature form with an integrand 
SJ(u, v)l sinM. As it happens, this parametrization is marginally 
but consistently more efficient than the Gaussian image, so com- 
bined with Lebedev quadrature, this is apparently the nominally 
fastest way to compute brightness integrals over ellipsoids. In 
practice, we found that the choice of surface parametrization has 
very little effect. 

4. Efficiency in lightcurve simulation and inversion 

We have shown above that, at an accuracy level of about 1% 
(corresponding to 0.01 mag in brightness), our scheme based on 
Lebedev quadrature is some five to ten times faster than triangu- 
lation, and orders of magnitude faster at higher accuracies. Thus, 
for direct brightness computations, i.e., lightcurve simulations, 
its efficiency is obvious. An important point is that functions are 
precomputed for chosen fixed points on S 2 just as with tessel- 
lation, so the procedure introduces no computational overhead. 
It is also easy to implement with Eqns. ( fTSl and ([T9T l using the 
known G(rf) or J(u, v) of the shape and any chosen scattering 
model S . The increase in speed is particularly important for the 
massive lightcurve simulations typically needed for estimating 
the performance of large-scale sky surveys such as LSST, Pan- 
STARRS or Gaia. The simulations can be carried out with shapes 
(i.e., G(t])) known for real asteroids from lightcurve inversion 
(Durech et al. 2010). Simulations with ellipsoids are fastest to 
run using Eqns. d28T i or (l22l . 

In lightcurve inversion, most of the computation time is spent 
on evaluating lightcurves and the corresponding partial deriva- 
tives dL/dPj for the shape and rotation parameters P,-, so the en- 
hanced speed in brightness computation yields an equally faster 
inversion procedure. The partial derivatives consume exactly the 
same computation time per evaluation point 77,- for both quadra- 
ture and tessellation. From Eq. ( [TBI we have 



^ WiG(T}i) 



dS(/i,iJo,a) 
dPs 



for rotation parameters Pr, and 



dL 
0P R 



dL sr 



dG( Vi ) 
dP R 



(29) 



(30) 



for shape parameters P$ ; the derivatives are of the same form 
as in the tessellation case ©. As before, functions (such as 
spherical harmonics) used in computing G(rf) are precomputed 
for chosen fixed points on S 2 just as with tessellation. Lebedev 
quadrature is thus easy to substitute for tessellation in the in- 
version procedure, and it introduces no computational overhead, 
so the difference in the computational cost between the two ap- 
proaches can be read directly from the A(N) (or A^(A))-curves in 
Fig. [5] We plan to include a Lebedev-based inversion version of 
the convexinv procedure downloadable in the software section 
of the DAMIT website (Durech et al. |20T0T >. 



The simplest way to utilize the superior speed of the quadra- 
ture approach in lightcurve inversion is the "black box" mode 
where the optimized G(rf) is used as an auxiliary unseen quantity, 
and the output parameters are those of rotation (pole direction 
and period). This is very useful in applications such as i) search- 
ing for the correct period in sparse photometric datasets from 
sky surveys when no period is evident a priori and a wide period 
range must be combed through (Kaasalainen 2004 , Durech et al. 
2006, 2007); ii) Monte-Carlo sampling to estimate the goodness- 
of-fit levels of the rotation parameters and thus their likelihood 
distributions; or iii) searching for and then fixing the pole and pe- 
riod prior to producing a high-resolution version of the shape as 
is usually done in the convex inversion procedure (Kaasalainen 
et al. 1200 11 1. Upon computing sample lightcurve inversion cases 
for targets previously analysed with tessellation, and using a re- 
duced number of Lebedev evaluation points in accordance with 
Fig. [5] we found that the shape results were virtually indistin- 
guishable from the previous ones, and the spin state parameters 
were essentially the same (e.g., typically within one or two de- 
grees of the previous pole solution), i.e., well within uncertainty 
limits and amounting to a slightly different initial guess. We em- 
phasize that the quadrature approach has very simple computa- 
tional plug-in properties as it has exactly the same form as, and 
is thus directly interchangable with, the tessellation part in the 
software, producing essentially the same results faster. 

Once the inversion procedure has provided the curvature 
function G (in, e.g., the form of the coefficients of a spherical 
harmonics series), this can be discretized into facet areas and 
normals usable by the Minkowski procedure by tessellation at 
the desired level of resolution, and this renders the result as 
a standard polyhedral model (Kaasalainen et al. |1992al [2001). 
Another way to use the Lebedev-based approach is to solve 
for the separate values G(rji) at Lebedev points 77, directly, in 
the same way that individual facet areas can be solved for 
in the tessellation-based inversion (using the exponential form 
G(r]i) - exp[a,] to guarantee positivity). Once these values are 
obtained, the facet areas and normals required by the Minkowski 
procedure can be defined by G(i]i)wi (when all w, are posi- 
tive) or by, e.g., using the Delaunay triangulation of the unit 
sphere shown in Fig. Q] The normals of the triangles give the 
facet normals of the final polyhedral model, and its facet areas 
are given by GAcr, where G is the average value of G at the 
three vertices of the corresponding Delaunay triangle, and Act is 
the area of the triangle. Thus the quadrature approach provides 
exactly the same low/high-resolution options as the standard 
tessellation-based procedure (e.g., convexinv) of Kaasalainen 
etal. (11992all200TT l. 

As described in Kaasalainen et al. (2001), the convex-shape 
procedure requires a regularization function (though usually 
with a very modest weight) to make sure that the resulting shape 
really is convex. This regularization is based on the three equa- 
tions fulfilled by convex surfaces: 



■In 



G(&,if/)Tn(8,il/)sm&d&dif/ = 0, 2 = 1,2,3; (31) 



i.e., a proper curvature function G has no first-degree terms Y"' 
in its representation as a spherical harmonics series. In regular- 
ization, the left-hand sides are enforced to be as close to zero as 
possible. The regularization integral is directly in the Lebedev 
form, so it is conveniently computed with the same quadrature 
as the lightcurve integral. 
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5. Conclusions and diSCUSSiOn Ostro, S., & Connelly, D. 1984, Icarus, 57, 443 

Wang X., & Carrington T. 2003, J. Theor. Comp. Chemistry, 2, 599 

This paper complements the study of the numerical computation 
of lightcurve integrals introduced and discussed in Kaasalainen 
et al. 01992a! 120011) . We have systematically investigated and 
compared different methods of computing lightcurve-type inte- 
grals over varying integration regions on the unit sphere S 2 . For 
smooth and continuous integrals over the whole of S 2 , Lebedev 
quadrature is usually the optimal method, but integrals over a 
part of S 2 depend on the case. For these integrals, the usual 
choice of (Gaussian) quadrature is computationally expensive 
with lightcurves because then i) either the integration limits are 
complicated or they can be transformed into simpler ones only 
at the price of complicating the computation of the integrand; 
and ii) when lightcurve computations are repeated for an ob- 
ject, functions have to be constantly re-evaluated, whereas for 
tessellation and for quadrature over the whole of S 2 they can 
be evaluated in one step. For lightcurve integrals, and other in- 
tegrals of similar nature, a modification of Lebedev quadrature 
by zeroeing a part of the quadrature sum outperforms both tes- 
sellation and quadrature over a part of S 2 (the latter when the 
required relative accuracy is not extraordinarily high - all real- 
istic applications fall into this category). The scattering function 
usually vanishes at the integration boundary, which increases the 
efficiency of the computation as it makes the nominal integrand 
continuous. 

Our scheme based on Lebedev quadrature is thus the opti- 
mal method for computating lightcurve integrals in practice. As 
a rule, Lebedev quadrature is five to ten times faster than trian- 
gulation at the accuracy level of 1% (0.01 mag), and it is orders 
of magnitude faster at accuracies of 0. 1 % and beyond (the slope 
of the accuracy curve is close to -2 rather than -1). The superior 
speed of the Lebedev approach is advantageous in any applica- 
tions in lightcurve inversion and simulation, and it is particularly 
suited to the automatic en masse analysis of thousands of tar- 
gets from large-scale surveys, such as Pan-STARRS or LSST, or 
for simulations requiring the computation of a large number of 
lightcurve points. 

Our method is very easy to substitute for tessellation in 
lightcurve simulation and inversion software since it is exactly 
similar to tessellation in its computational form: it merely uses 
fewer function evaluation points. The only thing that changes 
is that Lebedev weights (constants) are substituted for the fixed 
areas of facets on the unit sphere used in tessellation-based com- 
putations, and Lebedev points on the sphere are substituted for 
the fixed directions of the facet unit normals. 
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